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ABSTRACT 

In this paper we present the first formulation of the theory of non-linear particle acceleration in 
collisionless shocks in the presence of neutral hydrogen in the acceleration region. The dynamical 
reaction of the accelerated particles, the magnetic field amplification and the magnetic dynamical 
effects on the shock are also included. The main new aspect consists however in accounting for charge 
exchange and ionization of neutral hydrogen, which profoundly change the structure of the shock, as 
discussed in our previous work. This important dynamical effect of neutrals is mainly associated to 
the so-called neutral return flux, namely the return of hot neutrals from the downstream region to the 
upstream, where they deposit energy and momentum through charge exchange and ionization. We 
also present the self-consistent calculation of Balmer line emission from the shock region and discuss 
how to use measurements of the anomalous width of the different components of the Balmer line 
to infer the cosmic ray acceleration efficiency in supernova remnants showing Balmer emission: the 
broad Balmer line, which is due to charge exchange of hydrogen atoms with hot ions downstream of 
the shock, is shown to become narrower as a result of the energy drainage into cosmic rays, while 
the narrow Balmer line, due to charge exchange in the cosmic-ray-induced precursor, is shown to 
become broader. In addition to these two well-known components, the neutral return flux leads to the 
formation of a third component with intermediate width: this too contains information on ongoing 
processes at the shock. 

Subject headings: acceleration of particles - Cosmic Rays - Balmer emission - ISM: supernova rem- 
nants 



1. INTRODUCTION 

In the context of the so-called Supernova Remnant (SNR) paradigm for the origin of Galactic Cosmic Rays (CRs), 
the bulk of CRs are accelerated at SNR shocks through diffusive shock acceleration (DSA). Under this assumption, in 
order to explain the observed fluxes at Earth and, at the same time, the B/C ratio (sensitive to the grammage traversed 
by CRs during propagation) one is force d to require an accelera tion efficiency of 5 — 15% per SNR, depending on details 
of the propagation in the Galaxy (e.g.. iBlasi fe Amatdl2012aT fbl). In these conditions, the non- linear effects induced 
by the dynamical r eaction of accelerated p articles on the SNR shock cannot be neglected and reflect on the shape of 
the spectr um (e.g.. | Malkov &: Drurvll200lL for a review) and on the temperature of the ionized plasma downstream of 
the shock (|Decourchelle et al.ll2000D ~ Moreover, the streaming instability induced by accelerated particles ahead of the 
shock is thought to be res ponsible for the magnetic field amplification that is a crucial ingredient to reach maximum 
energies clos e to the kn ee (|Blasi. Amato fc"Ca prioli 2 0071 ) and to explain the thin non-thermal filaments observed in 
X-rays (e.g.,|VE3[2Qil). The non-linear theory of diffusive shock accelerat ion (NLDSA) including a ph e nomenological 
descr iption of all these effects has been presented in a series of papers (jAmato fc Blasil I2005L12006I : iCaprioli et all 
12009( 1. Both test-particle DSA and NLDSA lead to predict spectra which are rather hard at high energies, with a 
power-law index in mom entum < 4. This implies a steep dependence on energy of the Galactic diffusion coefficient 
(Berczhko & Volk 2007), if the CR spectrum at Earth is to be reproduced . Such energy dependen ce appears to be 
at odds with the observed anisotropy of Galactic CRs (jPtuskin et al.l [2006: Bl asi fc Amatol 120 1 2bl ). Accounting for 
the velocity of scattering centers might lead to somewhat steeper injection spectra thereby mitigating the problem 
(jCaDriolill20Ll . 

In general, anything that affects the structure of the shock also leads to changes in the spectrum of accelerated 
particles. One importan t source of shock m odification is represented by the presence of neutral atoms in the acceleration 
region: as discussed by IBlasi et al.l ([20I2I . hereafter Paper I), the reactions of charge exchange (CE) and ionization 
produce important effects in the shock region. The most important of these effects is induced by the so-called neutral 
return flux (NRF): a hydrogen atom crossing a collisionless shock front can charge exchange with a hot ion downstream 
of the shock, thereby leading to the formation of a cold ion and a hot neutral. It may happen that the latter is produced 
with velocity in the direction of the shock; being neutral, it can recross the shock and suffer CE or ionization upstream, 
hence releasing energy and momentum there. This important phenomenon leads to heating of the upstream ionized 
plasma and to a decrease of the Mach number. Even in the context of test-particle DSA it is easy to show that 
this phenomenon leads to a steepening of the spectrum of particles accelerated in a SNR where appreciable neutral 
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hydrogen is present (Paper I). The NRF affects the upstream plasma on a scale of a few interaction lengths of CE. 

The presence of neutral hydrogen also represents a unique diagnostic tool of the acceler ation proces s, through the 
Balmer emission associated with CE to excited states: in a Balmer dominated shock (sec lHen3 [20091 for a review 
the Balmer line has a narrow and a broad component, the first associated with hydrogen atoms that did not suffer 
CE with high er temperature ions, and the second that reflects the te mperature of shocked ions, due to CE reactions 
downstream (j Chevalier fc Raymond! [T9781 Ivan Adelsberg ct al J 119781 ) . In the presence of CR acceleration the ions' 
temperature downstream is lower, and this should reflect in decreased Full- Width at Half-Maximum (FWHM) of 
the broad Balmer line. If the CR pressure is large enough, a CR-induced precursor is formed in front of the shock 
(jMalkov fc Drury|[200l so that CE can also take place upstre am and eventually lead to a broadening of the narrow 
Balmer line. In addition, it was shown in lMorlino et al.l (|2012t hereafter Paper II), that the shape of the Balmer line 
is also affected by the NRF discussed above through the formation of an intermediate component of the line. 

Remarkabl y, anomalous width s of both broad and narrow lines have been observed in Balmer-dominated shocks. 
For instance. iHelder et al.l (|2009[ ) combined proper-motion measurements of the shock and broad Ha line width for the 
remnant RCW 86 to demonstrate that the inferred low temperature behind the shock may be interpret ed as a signal 
that a sizable fraction of the energy is being channeled into CRs. A similar conclusion has been reached bv lHelder et al.l 
(pOlOf ) by analizing the SNR 0509-67.5 in the LMC, even if in this case the shock velocity has been estimated using 
an evolution model with information derived from the line width of elements in the shocked ejecta (wich can be used 
to infer the speed of reverse shock). 

Anomalous width of narrow Balmer lines have been also observed in several SNRs (see, e.g. iSollerman et al.ll2003() . 
The width of such lines is in the 30-50 km s _1 range, implying a pre-shock temperature around 25,000-50,000 K. If 
this were the ISM equilibrium temperature there would be no atomic hydrogen, implying that the pre-shock hydrogen 
is heated by some form of shock precursor in a region that is sufficiently thin so as to make collisional ionization 
equilibrium before the shock unfeasible. The CR precursor is the most probable ca ndidate to e xplain such a broadening 
of the narrow line width, though other mechanisms have also been proposed fsee lHenell2009l . for a revi ew) . 

A p revious attempt to include neutral particles in the shock acceleration theory has been done by Morlino et al.l 
using a fluid approach. Nevertheless, as discussed in Paper I, neutrals cannot be treated as a fluid be cause 
the length scales involved are much smaller than the equilibration scale. More recently iRavmond et al.l (|2011[ ) tried 
to describe the interaction of neutrals with the CR precursor using a phenomenological approach where both the CR 
spectrum and spatial distribution are assumed rather than calculated. 

An accurate description of the CR acceleration process in the presence of neutrals can only be achieved by formulating 
a NLDSA theory where the dynamical action of neutrals is taken into account. This is a needed preliminary step, if 
one wants to make use of measurements of Balmer line emission from SNRs to quantitatively assess their efficiency 
as particle accelerators. In this paper we present the first treatment of NLDSA in partially ionized media and discuss 
how to use the theory to calculate the shape of the Balmer line and infer information on particle acceleration at the 
shock. 

The paper is structured as follows: in Section [2] we describe the calculations used to describe neutrals, ions and the 
acceleration process and we explain in detail how to assemble the parts together to describe NLDSA in the presence of 
neutrals. The formalism previo usly introduced in Paper I — as far as neutrals are concerned — and in (jAmato fc Blasil 
120051 [20061 : iCaprioli et al.ll2009[ ) — as far as NLDSA is concerned — will be heavily used. In Section [3] we illustrate our 
results in terms of shock modification, spectrum of accelerated particles and shape of the Balmer line. We summarize 
and conclude in Section @] 

2. CALCULATIONS 

We consider a plane-parallel shock wave propagating in a partially ionized proton-electron plasma with velocity 
V^h along the z direction. The fraction of neutral hydrogen is fixed at upstream infinity where ions and neutrals are 
assumed to be in thermal equilibrium with each other. The shock structure is determined by the interaction of CRs 
and neutrals with the background plasma. Both CRs and neutrals profoundly change the shock structure, especially 
upstream where both create a precursor: the CR-induced precursor reflects the diffusion properties of accelerated 
particles and has a typical spatial scale of the order of the diffusion length of the highest energy particles. The neutral- 
induced precursor develops on a spatial scale comparable with a few interaction lengths of the dominant process 
between CE and ionization. The downstream region is also affected by the presence of both CRs and neutrals and the 
velocity gradients that arise from ionization have a direct influence on the spectrum of accelerated particles. 

A self consistent description of the overall system requires to consider four mutually interacting components: thermal 
particles (protons and electrons), neutrals (hydrogen), accelerated protons (CRs) and turbulent magnetic field. We 
neglect the presence of helium and heavier chemical elements. In the following subsections we write down the basic 
equations used to describe each component, including the interaction terms. We are looking for stationary solutions, 
hence all equations arc written as time-independent. The interaction terms make the system highly non lin ear and the 
solution is found usin g a iterative scheme similar to the one we introduced in some of our previous work (IBlasil 120021: 
lAmato fc~ Blasi 2005], Paper I). The details on how to generalize such scheme to the case in which neutrals are present 
are discussed in § 12.51 

2.1. Neutral distribution 

The behavior of neutrals has been extensively discussed in Paper I. Neutral hydrogen interacts with protons through 
CE and ionization and with electrons through ionization only. The hydrogen distribution function, fN(v,z), can be 
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described using the stationary Boltzmann equation 

df N (v,z) 



p N f i (v,z)-[p i + p e }f N (v,z), (1) 



where z is the distance from the shock (which is located at the origin), v z is the velocity component along the z 
axis and the electron and proton distribution functions, fi(v,z) and f e (v,z), are assumed to be Maxwellian at each 
position. The collision terms, Pkfi, describe the interaction (due to CE and/or ionization) between the species k and 
/. The interaction rate Pk is formally written as 

p k (v,z) = I d 3 wv re ia(v re i)f k (w,z) , (2) 



where v re i = \v — w\ and a is the cross section for the relevant interaction process. More precisely, Pn is the rate of 
CE of an ion that becomes a neutral, Pi is the rate of CE plus ionization of a neutral due to collisions with protons, 
while (3 e is the ionization rate of neutrals due to collisions with electrons. A full description of the cross sections used 
in the calculations can be found in Paper II. 

A novel method to solve Eq. ([T]) has been developed in Paper I, where we showed how to calculate /at starting from 
the distribution of protons and electrons. The method is based on decomposing the distribution function /jv in the 
sum of the neutrals that have suffered 0, 1, 2, ... , k processes of charge exchange. Each distribution function is named 

fx , and clearly /jv = J2T=o /jv • Using this formalism, Eq. ([T]) can be rewritten as a set of k + 1 equations, one for 
each component: 

v z d z f ( x ] = - \fii + p e ] ffi forfc = 0, and (3) 
vAfP = p ( t l) fi - [pi + p e ] forfc>0. (4) 

(k) . (k—1) 

Each one of these equations is linear, and the solution for can be obtained from . The boundary conditions 

are imposed at upstream infinity where ions and neutrals are assumed to start with the same bulk velocity and 
temperature, therefore charge exchange occurs at equilibrium and the distributions do not change. As a consequence, 

the boundary conditions are fffl (v, — oo) = (njv,oAHo)/i(v, — oo) and /Jy (v, — oo) = 0. A good approximation to the 

total distribution function of neutrals is obtained when a sufficient number of is taken into account. This number 
is determined by the physical scale of the problem. A rough estimate can be obtained comparing the lengthscale for 
CE whith the maximum between the ionization lengthscale and the precursor length. In general, for reasonable values 
of the ambient parameters, the longest scale, and the one relevant for comparison, is that of the CR-induced precursor. 

2.2. CR distribution 

The isotropic distribution function of CRs satisfies the following transport equation in the reference frame of the 
shock: 

d_ 

dz 



D(z,p)-^ 



df c ldu df c 

U ^ + 3d-z P W + Q ^-°- (5) 



The z-axis is oriented from upstream infinity (z = —oo) to downstream infinity (z — +oo) with the shock located at 
z = 0. We assume that the injection occurs only at the shock position and is monoenergetic at p = p ln j, hence the 
injection term can be written as Q(z,p) = Qo(p)S(z) where Qo(p) = (rji n jni/4Trpf n j)S(p — pi n j). Here m is the number 
density of ions immediately upstre am of the subshock , and ?7i n j is the fraction of particles that is going to take part in 
the acceleration process. Following iBlasi et al.l (|2005l ). ?7inj can be related to the subshock compression factor as: 

r/inj = 4/ (3 v^F) (i? sub - 1) efnje"^ ■ (6) 

Here £i n j is defined by the relation p; n j = £i n jPt/i.2 , where Pth,2 is the momentum of the thermal particles downstream. 
£i n j parametrizes the poorly known microphysics of the injection process and is taken as a free parameter with a typical 
value between 2 and 4. 

The diffusion properties of particles are described by the diffusion coefficient D(z,p). We assume Bohm diffusion in 
the local amplified magnetic field: 

D(z,p) = ±cr L [6B(z)}, (7) 
where r^{8B) = pc/[eSB(z)] is the Larmor radius in the amplified magnetic field. The calculation of SB is described 



The solut ion o f Eq. ([51) can be easily obtained by using standard techniques. Here we generalize the approach of 
iBlasil (|2002l ) and lAmato fc Blasil (|2005l ) in order to relax the assumption of null gradient of the distribution function 
downstream (homogeneity), which does not apply to the general case in which CRs and neutrals are present. 
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The solution of Eq. ([5]) in the upstream is the s ame found by 



approximated by the following expression (see also iCaprioli et al 



Amato fc Blasil (120051 ) . Such a solution can be well 
120 101) : 



fu P {z,p) = / (p)exp 



D(z',p) 



(8) 



where fn(jp) = f (0, p ) is the distribution function at the shock position. Through the same approach used by 
lAmato fe Blasil ()2005l ) we can express the solution in the downstream in the following implicit form: 



, , x f ( x , 1 f Z dz' f°° , „du Of f 

w*,p) = / (p) + 3 y o ^ d - ^^ ex p [-j z 

We found that an excellent approximation to the expression in Eq. ^ is given by: 
fdown{z,p) = / (p)exp 



15 dy 



| Qo(p) 









az cxp 



-dz" 



where qo(p) = — dhi fo(p) /dlnp is the slope of the distribution function at the shock position. Both Equations 
(fTU|) give f(z,p) = fo(p) when the downstream speed it(z) is constant. 
The distribution function at the shock is obtained, as usual, by integrating Eq. ([5]) across the subshock: 



Ux - u 2 df 
-^- P ~dp- 









az 


i 


az 



+ Qo(p) 



(9) 

(10) 
and 

(11) 



Here and in the rest of the paper, labels "1" and "2" refer to quantities immediately upstream and immediately 
downstream of the subshock, respectively. The value of [Ddf/dz] 1 can be obtained by integrating Eq. ([5]) between 
upstream infinity (z = — oo) and position 1 (|Blasi 2002). In the standard stationary solution the downstream distribu- 
tion is homogeneous, hence [D df /dz] 2 = 0. This condition docs not hold in our case because of ionization of neutrals 
downstream. The general expression can be derived from Eq. (j9]) and reads 



D 



dz 



-kip) 



qoip) 



dz 



du 



-dz" 



Substituting the last expression in Eq. (1111) we obtain for fo(p): 



kip) 



r)n 



3i?t 



47Tpf nj RtotUpip) - 1 



exp ■ 



3R 



tot 



Pinj p' RtotU p (p') 



D 



(12) 



(13) 



where we defined U p (p) = u p /V s h and u p {p) can be interpreted as an estimate of t he mean plasma velo city that a 
particle with momentum p experiences in the precursor region (see Equation (8) in lAmato fc Blasi 2005). Eq. (jlc 
differs from the standard solution only for the term including K(p), where 



K(p) = —— I dz—exp 



dz 



^dz< 
D 



(14) 



The mean plasma speed experienced in the downstream by particles that are able to return to the subshock is 
U2 + V s hK(p). When the do wnstream plasma has a constant velocity, K{p) — » and Eq. (flU)) reduces to the standard 
solution (sec Equation (7) in lAmato fc B lasi 2005|). 

2.3. Transport equation for waves 

Following I Amato fc Blasil (|2006| ) we describe the scattering of accelerated particles in the acceleration region as due 
to the waves that the particles generate through resonant streaming instability. These waves are also damped due to 
several processes. In particular, when the plasma is not fully ionized, the presence of neutrals can damp Alfvcn waves 
via ion-neutral damping. 

The equation for transport of waves can be written as: 



8F w (k,z) . ,dP w (k,z) 



dz 



P w (k,z) [o- C Rik, z) 



TH 



(k,z)] 



(15) 



where F w and P w are, respectively, the energy flux and the pressure per unit logarithmic bandwidth of waves with 
wavenumber k. a is the growth rate of magne tic turbule nce, while Tth is the damping rate. For resonant wave 
amplification the growth rate of Alfven waves is (IBelllll978D : 



cr C R(fc,x) 



47T va(x) 
3 Pw i^i %) 



p 4 v(p) 



dx 



(16) 



p=p(k) 
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where p = p(k) = eB/km p c is the resonant momentum. The damping of the waves is mainly due to non-linear Landau 
damping and ion-neutral damping. For the sake of simplicity here we adopt a phenomcnological approach in which the 
damping results in a generic Alfvcn heating at a rate Tth = '(th^cr- This expression assumes th at a fraction ?)th of 
the p ower in amplified waves is locally damped and results in heating of the background plasma ( Berezhk o fc Ellison! 
[1993). 

Following ICaprioli et all (pOOl . Eq. <^ can be solved by integrating in /c-space and using the relation between 
energy density and pressure of Alfven waves, namely F w = 3P w u, which holds upstream when va *C V s h- This leads 
to the following equation: 



W{X) &M = Va {z) [1 - , T „] d -^l - ZP w ( Z ) dU{z) 



(17) 



dz dz dz 

u/V sh , V A = v a /V Bh , P = P/(poV^ h ). Eq. dTTJ) can be used to derive P, 



where we used the normalized quantities U 
in the upstream region. 
The jump conditions for the magnetic field are calculated following ICaprioli et al.l (|2009j) 



F 

I! 



(18) 
(19) 



In the downstream section we solve the same transport Eq. (|15|) but we assume the relation between energy flux and 
pressure that derives from Eqs. (fl~8| and ([T9]) . namely F w = (1 + 2R su \>)uP w . Moreover, in the downstream region we 
assume for simplicity that damping occurs on scales much larger than the typical scale length for CE and ionization 
(this assumption is not strictly needed, it simply makes the results of the calculations easier to interpret). 

2.4. Dynamics of the ionized gas 

The dynamics of the background plasma is affected by the presence of accelerated particles and by CE and ionization 
of neutrals. Protons and electrons in the plasma are assumed to share the same local density, pi{z) = p e (z), but not 
necessarily the same temperature, i.e., Ti(z) may be different from T e (z). The equations describing the conservation 
of mass, momentum and energy taking into account the interactions of the plasma fluid with CRs are: 



9 I 

■X- [PiUi 

oz 



+ p N ] = , 



dz 



d_ 

dz 



2 PiUi 



PtU, 



P„+Pr + P, 



Pn] 



IgPgUi 



7s 



1 



F, 



N 



= o, 

dPc 
dz 



(20) 
(21) 

rPu, . (22) 

■f w^)/jv are respectively the fluxes of 



Here p^ = m# J d 3 vv^f^, P^ = mu J d?"vv^f^ and F^ = m#/2 J d 3 vv^(v^ 

mass, momentum and energy of neutrals along the z direction (labelled as ||). They can be easily computed once the 
neutral distribution function is known. P w and F w are the pressure and energy flux of waves as defined in section 
§ 12.31 while P c is the CR pressure computed from the CR distribution function: 



Pc{z) 



4tt 



dpp 3 v{p)f c (z,p) . 



(23) 



The dynamical role of electrons in the conservation equations is usually neglected due to their small mass. However, 
collective plasma processes could contribute to equilibrate electron and proton temperatures, at least partially. If the 
equilibration occurs in a very efficient manner, the electron pressure cannot be neglected and the total gas pressure 
needs to include both the proton and electron contributions, namely P g = Pi + P e = Pi(\ + j3), where j3(z) = T e /Ti 
is the electron to proton temperature ratio and is taken here as a free parameter. We note that the electron-ion 
equilibration level can be different between upstream and downstream because the plasma conditions in these two 
regions are different. Hence we distinguish between upstream and downstream using two separate parameters, /3 up and 
/?down, respectively. On the right side of Eq. (|22| we include the heating produced by the compression of CR acting 
on the plasma plus the heating due to the damping of magnetic waves. 

The computation strategy we adopt is as follows: we first solve Eqs. (|20p - d22l) for the upstream dynamics starting 
from the boundary conditions at upstream infinity, i.e. Pi(— oo) = p^o and Ui(— oo) = V^h- In order to determine 
the conditions immediately downstream of the subshock we use the jump conditions derived from the conservation 
equations: 



[piUi + fi N } 1 = , 

[P^ + Pg + P W }\ =0, 



lg p g u i 
7»-l 



F, 



0. 



(24) 
(25) 



(26) 
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Fig. 1. — Flow chart of the iterative method used to find the solution. 



where the brackets indicate the difference between quantities upstream and downstream of the shock ([X]f = X2 — X1). 
The terms describing neutrals and CRs are absent because of the continuity of Jn and f c across the subshock. Finally 
we solve the conservation equations, Eqs. ([2"Ul) - (f2"2")) in the downstream region, using the values of the quantities at 
position 2 as boundary conditions. 

2.5. Iterative method for the solution of the non-linear problem 

In order to solve the set of non-linear equations involving neutrals, ions, CRs and magnetic field, we adopt an 
iterative method that we illustrate in Fig. (|T|) in the form of a flow diagram. The input quantities are the values of 
the shock velocity and all environmental quantities at upstream infinity, where the distribution function of neutrals 
is assumed to be Maxwellian at the same temperature as that of ions. The method proceeds through two nested 
cycles of iteration: in the first cycle, we start with a distribution function of ions that corresponds to an unperturbed 
distribution that one would obtain in the absence of both neutrals and CRs and we calculate the neutral distribution 
function /jv(z) from Eq. ([1]). With the given /jv(^) we iterate to obtain the ion distribution, the CR distribution 
and the wave energy density. The outer loop of iteration allows us to update the neutral distribution function using 
the newly available ion and CR distribution functions. The iteration ends when the predefined accuracy is achieved. 
At this point the distribution functions are used to calculate the structure of the Balmer line, following the method 
illustrated in Paper II. 

The typical computation time required to complete the calculations for a given set of environment parameters is of 
the order of one hour on a standard laptop. The most time-consuming part is the solution of Eq. ([lj which increases 

linearly with the number of partial functions, /J^ , required for a correct description of the neutral distribution function. 

(k) 

In particular the number of increases linearly with the ratio between the precursor length and the CE interaction 
length. 

3. RESULTS 

The formalism presented above represents the first theory of particle acceleration at collisionless shocks in the 
presence of neutral atoms and allows us to calculate the shock structure, the spectrum of accelerated particles and the 
Balmer emission from the shock region. In the subsections below we describe the main results concerning these three 
topics. 

3.1. Shock structure and spectra of accelerated particles 

We adopt a test case consisting of a shock moving at V s h = 4, 000 km s _1 in a medium with temperature T = 10 4 K, 
magnetic field B = 10 fiG, and densities p^o = Pn.o = 0.05 cm -3 (the neutral fraction is 50%). CRs are injected with 
an efficiency of ~ 20% (in fact we fix the injection parameter £ = 3.7) and the efficiency of Turbulent Heating (TH 
hereafter) is fixed as 7?th = 0.5 (see Eq. [T7|) . The shock structure is illustrated in Fig. [2j the top panel shows the 
pressure of gas (continuous solid line), magnetic field (dotted line) and CRs (dashed line) as functions of z, plotted 
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Fig. 2. — Example of spatial profiles of relevant quantities for a shock propagating in partially neutral plasma. Top panel: pressure of 
ions, CRs and magnetic field. Middle panel: ion temperature. Bottom panel: ion velocity. Quantities are shown in the shock frame with 
the subshock located at the origin. The z axis is shown in logarytmic scale both upstream and downstream and zj^ax = £>(pmax)/V^h is the 
diffusion length of particles with the maximum momentum. These profiles are computed for: V s h = 4, OOOkms - 1 , To = 10 4 K, Bo = 10 fiG, 
PiO = PN,0 = 0.05cm -3 . The injection parameter for CRs is £i n j = 3.7, which corresponds to an acceleration efficiency ~ 20%. 



here in logarithmic scale for both the upstream and downstream regions. All pressure components are normalized 
to the ram pressure Pi^V^ at upstream infinity. Clearly the CR pressure is continuous across the shock front. The 
self-generated magnetic field acquires a pressure that is of the same order of magnitude of the therma l pressure, which 
implies that its dynamical effect is non- negligible, as already discussed in lCaprioli et al.l (|2008l [2009f ). 

The background ions are substantially heated in the precursor due to TH, as it is shown by the increase in the gas 
temperature on a spatial scale comparable with z* lax = z*(p max ) = D(p max )/V a h, i- e -> the diffusion length of particles 
at Pmax (mid panel of Fig. [2]) . The upstream profile of the fluid velocity, showed in the bottom panel of Fig. [2] as 
normalized to V s h, is a distinctive signature of the CR-induced precursor. 

Finally, one last thing to notice in this Figure is the increase of both the gas pressure and velocity at some distance 
from the shock in the downstream. This is a result of progressive ionization of fast cold neutrals, and indeed occurs 
on scales which are of the order of the ionization length. We will see that the situation is different, with the plasma 
velocity decreasing, rather than increasing, towards downstream infinity, when the shock velocity is lower and the 
NRF starts to play an important role (see Fig. |4] below) . Spatial variations of the plasma velocity are important for 
determining the accelerated particle spectrum at the shock, whose slope at each given energy always reflects the mean 
compression ratio experienced by particles of that energy while moving away from the shock. 

The spectrum of accelerated particles at different locations upstream and downstream of the shock is shown in Fig. [3] 
for a neutral fraction = 0.5. The solid lines refer to the downstream section of the plasma, the thick dot-dashed 
line is the spectrum at the subshock location, the dashed lines are the spectra at different locations upstream of the 
shock. 

One can see that, moving away from the shock in the upstream direction, the spectrum is more and more depleted of 
the lowest energy particles: at any position z, only particles energetic enough to have a diffusion length D(p)/V s h ~ z 
are present. The transport of low-energy particles in the downstream, is rather dominated by advection, and the 
spectrum is never truncated: at downstream infinity, the low-energy spectrum has the same slope as at the shock, but 
with a higher normalization. This is due to the compression that occurs around 0.1 < z/z* lax < 1, where ionization 
takes place. At the highest energies, instead, diffusion dominates over advection and the CR spectrum far downstream 
is expected (and found) to be exactly the same as at the shock location. As a result, at large distances downstream of 
the shock the spectrum turns out to be slightly steeper than at the shock. 

As already illustrated in Paper I in the test-particle limit, the effect of the NRF is to create an upstream precursor 
that steepens the spectra of accelerated particles. For this reason the spectra of accelerated particles are less concave 
than they would be in the presence of the CR reaction alone. This result is qualitatively general, though the strength 
of the effect is quantitatively dependent upon the specific environment that the shock propagates into, as discussed 
below. 

In Fig. |4] we show the normalized plasma velocity as a function of the distance from the shock for different values 
of the neutral fraction (left panel), and the slope of the spectrum of accelerated particles as a function of momentum 
(right panel). All the calculations refer to a shock speed V s h = 2, OOOkms" 1 , a total gas density at upstream infinity 
no = 0.1 cm~ 3 , an unperturbed magnetic field strength Bq = 10 fiG and an injection parameter £i n j = 3.7. The solid 
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Fig. 3. — CR momentum spectrum multiplied by p 4 for different distances from the sub-shock: the dash-dotted (black) line is the 
spectrum at the sub-shock position, the dashed (red) lines show the spectrum in the upstream and solid (blue) lines show the spectrum 
downstream. The position z is indicated for each curve in units of z*^^. 
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Fig. 4. — Left panel: Velocity profile of the upstream and downstream fluids normalized to the shock velocity, as a function of the distance 
from the shock. The shock is located on the left side and the (absolute) distance is in unit of z*(p m ax). The profiles arc shown for a given 
shock speed V^h = 2, 000km s _1 and for different neutral fractions /ijv from up to 80% as in the legend. Other parameters are £j n j = 3.7, 
Bq = 10 fiG and no = 0.1 cm -3 . Right panel: Slope of the spectrum of accelerated particles for the same shock configurations as in the left 
panel. 



lines refer to the case with fully ionized plasma (hp, = 0). In this case the shock is appreciably modified by the CR 
pressure, the acceleration efficiency being ~ 30%, and the spectrum is quite concave. The spectrum is steeper than 
p~ 4 only for energies below ~ 10 GeV. 

Due to the relatively-low shock velocity, the effects of CE and ionization are quite relevant; therefore, we expect the 
results to change appreciably for non vanishing values of h^, due to the NRF. This is illustrated by the long-dashed, 
dashed and dotted lines in Fig. 0] For hp, = 0.2 (long-dashed line) and hp, = 0.5 (dashed line) the formation of 
the precursor upstream is dominated by the dynamical reaction of CRs on the plasma, although the heating of the 
upstream due to the NRF is very evident. In these two cases the acceleration efficiency is ~ 20% (for hp, = 0.2) and 
~ 15% (for hp, = 0.5). The spectra in the same cases are correspondingly less modified by the CR pressure (less 
concave) and, on average, somewhat steeper: for hp, = 0.2 the spectrum is steeper than p -4 for energies < 100 GeV, 
while for h p, = 0.5 it remains steeper than p~ 4 even at E ~ 10 4 GeV. 

The case with large neutral fraction, hp, = 0.8 is rather extreme: the precursor becomes very short, as it becomes 
dominated by the NRF, and the acceleration efficiency is still ~ 15%. The spectrum is quite steeper than p~ 4 at 
all energies, so that the bulk of the pressure is carried by relatively low energy particles. This implies that the CR 
pressure can slow down the incoming plasma only very close to the subshock (see left panel in Fig. 3]) . 

As discussed in Paper I, the effects of the NRF are most important for shock speeds V s h < 3, OOOkms -1 ; for larger 
shock speeds, the typical relative velocity between neutrals and ions is such that the cross section for CE is so small 
that atoms get ionized before suffering a CE reaction. As a consequence, the NRF is suppressed. It is useful to explore 
the dependence of the spectrum on the shock velocity in a case in which the neutral fraction is h p, = 0.5. In the 
left panel Fig. [5] we plot the normalized plasma velocity as a function of distance from the shock for shock velocities 
V sh = 2,200kms" 1 (solid line), V sh = 3, OOOkms" 1 (dashed line) and V^h = 4,000kms- 1 (dotted line). The other 
parameters are the same as in Fig. |4] The right panel shows the spectral slope for the same values of the shock velocity. 
For comparison, the thick lines refer to = 0.5, while the thin lines of the same type refer to completely ionized gas 
(hi, = 0). 

The left panel qualitatively confirms our expectations based on the test-particle results: when increasing the shock 
speed from 2, 200 kms -1 to 4, 000 kms -1 , the effect of neutrals becomes smaller, namely the shock becomes increasingly 



Collisionlcss shock in a partially ionized medium: 



III. Efficient cosmic ray acceleration 



9 




more modified by the reaction of CRs. This is even more evident in the right panel of Fig. [3] for V a h = 2, 200 km s" 1 
a neutral fraction of 50% is sufficient to make the spectrum of accelerated particles steeper than p~ 4 , at all energies 
(thick solid line). With the same shock speed, but for a completely ionized plasma, the spectrum is instead harder 
than p~ A above ~ lOGeV (/ijv = 0, thin solid line). 

3.2. Balmer emission 

The Balmer line is affected by the presence of CRs in at least two ways. First of all, the width of the broad Balmer 
line is expected to be narrower if efficient CR acceleration takes place. This is a straightforward consequence of energy 
conservation: energy is channelled into accelerated particles, therefore the plasma temperature downstream is lower. 
Moreover, the formation of a CR-induced precursor leads to a non-negligible heating upstream of the subshock, so that 
the narrow Balmer line may become broader. In this section we investigate these possibilities in a more quantitative 
way, by using the formalism introduced in Paper II. The different components of the line are identified as discussed in 
Paper II. The global shape of the Balmer line is fitted with three Gaussian components: a narrow line, an intermediate 
line (reflecting directly the effect of the NRF) and a broad line. 

In Fig. [5] we show the temperature throughout the shock region for different values of the efficiency of TH, ?7th, for 
a shock propagating at 4, 000 km s" 1 in a medium with density 0.1 cm -3 . The case without CR acceleration is shown 
as a thick solid (black) line. In this case the precursor is solely due to the NRF and is limited to the region very close 
to the subshock. 

When the acceleration process is turned on, several phenomena become visible. In the absence of TH (t/th = 0, 
thin solid line), the upstream fluid is mildly heated by the adiabatic compression due to the CR-induced precursor 
only. Nevertheless, the temperature of the downstream plasma drops by almost a factor 2 as a consequence of the 
energy channelled into the accelerated particles. Increasing the value of t?th, the heating in the precursor becomes 
more and more marked, and it extends over larger and larger distances. However, the downstream temperature is not 
appreciably affected by a change in the efficiency of TH. 

The shape of the global Balmer line emission is plotted in the left panel of Fig. while the right panel shows a 
zoom in the region of the narrow Balmer line. The curves refer to the same cases as in Fig. [31 labelled as indicated. 
The width of the broad component of the Balmer line is appreciably reduced when CR acceleration is turned on, but 
the results are not sensitive to the amount of TH. On the other hand, as one can notice in the right panel, the narrow 
Balmer line becomes broader when CR acceleration is efficient; the effect is more pronounced when non-negligible TH 
is taken into account. Hence, the distribution function of neutrals becomes broader mainly because of the scattering 
with a warmer ion distribution in the far precursor (i.e., where TH is more effective), rather than because of the NRF, 
which operates only within a few CE interaction lengths from the subshock. 

In Fig. [5] we show the FWHM of the broad Balmer line as a function of the CR acceleration efficiency ecR (left 
panel) and as a function of the shock velocity V s h (right panel). In the left panel, the two curves illustrate the 
change in the FWHM obtained by assuming full or partial temperature equilibration between electrons and protons 
in the downstream (/?down = 1 and /3down = 0.01 , with /3down the ratio between electron and ion temperature), for 
V sh = 4, 000 km s -1 and n = 0.1 cm~ 3 . 

The FWHM decreases when the energy density in CRs increases, confirming the very important notion that the 
width of the broad Balmer line can be used to measure the CR acceleration efficiency. There is however an important 
caveat to keep in mind: for a given value of ecR; the uncertainty in the electron-proton equilibration translates into 
a spread in the FWHM of ~ 500 km s" 1 . Hence, unequivocal evidence of efficient CR acceleration is likely to be 
achievable only if the measured FWHM is below the line corresponding to full equilibration. Even in this case the 
value of ecR is likely to be quite uncertain because of the unknown level of equilibration. For instance, if a FWHM of 
2, 200 km s -1 were measured for a shock with V^h = 4, 000 km s -1 , the result could be interpreted as a measurement of 




Fig. 6. — Temperature of the plasma in the shock region for V s h = 4, 000 kms 1 , no = 0.1 cm 3 , hjv = 0.5 and Bo = 10 fj,G. When CR 
acceleration is turned on, we assume £inj = 3.5 and p max = 50 TeV/c. The thick solid (black) line refers to the case with no CRs. The thin 
(red) solid line refers to the case with CRs but no TH. Different values of «th are as indicated. 




Fig. 7. — Volume-integrated line profile of Balmcr emission for the same cases shown in Fig. [6] The right panel shows a zoom on the 
narrow line region, i.e. the shadow box in the left panel. The thick solid line shows the case without CRs, while different thin lines are 
calculated at fixed £j n j = 3.5 for different values of the TH efficiency. 



£cr ~ 0.25 in the case of full equilibration. In the more general case of partial equilibration, the above value of ecu can 
be taken as a lower limit. On the other hand, for the same shock velocity, if one measured a FWHM of 2,600kms~ 1 , 
it would be hardly possible to say anything at all about the efficiency of CR acceleration, without knowing the level of 
electron-proton thcrmalization. In order to have an independent estimate of /3, it would be very helpful to pair optical 
observations with the emission of SNR shocks in the (thermal) X-rays as well. 

In the right panel of Fig. [8] we fix the injection parameter £; n j = 3.5 and we calculate the FWHM of the broad 
Balmer line as a function of the shock velocity. The two thin (red) lines refer to the cases with no CR acceleration with 
/?down = 0.01 (upper curve) and /3down = 1 (upper curve). The dashed (t/th = 0) and dash-dotted lines (tjth = 1) refer 
to the cases with CR acceleration (efficiencies in parentheses), again for /?down = 0.01 (upper curves) and ftown = 1 
(lower curves). From the plot it is clear that the TH does not affect appreciably the width of the broad line. As pointed 
out above, if the electron-proton thermalization is inefficient, it may be difficult to draw any strong conclusion about 
the CR acceleration efficiency unless additional information is available to constrain the value of /3down- Nevertheless, 
measuring a FWHM of the broad line below the solid thick line allows one to put a lower limit on the CR acceleration 
efficiency. 

Additional discrimination power might derive from a combined analysis of the broad and narrow Balmer line com- 
ponents. The FWHM of the narrow Balmer line as a function of the maximum momentum p max is plotted in Fig. [5] 
for 77TH = 0.2 (top left panel), ryxH = 0.5 (top right panel) and 77th = 0.8 (bottom panel). The three curves in each 
panel are obtained for £j n j = 3.5, 3.7, 3.8. The shock velocity and total density are set to V s h = 4, 000 kms" 1 and 
no = 0.1 cm -3 , which approximately correspond to 6cr = 0.4, 0.2 and 0.1 (weakly dependent upon ?7th) for the three 
values of £i n j. The maximum momentum, p max determines the spatial extent of the CR-induced precursor. Larger 
values of p max imply that there is more time (space) for depositing heat in the upstream, and the width of the narrow 




Fig. 8. — FWHM of the broad Balmer line as a function of the CR acceleration efficiency ecR {left panel) and as a function of the shock 
velocity V s h (right panel). The upper set of curves (thin lines) correspond to the case ftiown = 0.01, while the lower curves (thick lines) 
correspond to ftjown = !• ln the right panel, the values near the points refer to the associated CR acceleration efficiency. 




Fig. 9. — FWHM of the narrow line as a function of the maximum momentum of accelerated protons. The three panels refer to ?7th =0.2, 
0.5 and 0.8. Each panel shows three lines calculated assuming different injection parameter, ginj =3.5, 3.7 and 3.8 which approximately 
correspond to ecu = 0.4, 0.2 and 0.1. 



Balmer line broadens correspondingly. The effect becomes more pronounced for larger values of the parameter ryxH- 
It is worth noticing that the behavior of the curves in Fig. [9] is the same in all cases and reflects the slight broadening 
of the narrow Balmer line as due to the NRF. In fact, one can estimate the momentum at which the diffusion length 
equals the size of the precursor induced by the NRF, which is of order Anrf ~ 10 16 cm for the parameters used here. 
By imposing D(p max ) /V a h = PmaxC 2 /3e SB = Anrf one gets p max ^5—10 TeV/c (from our calculations the amplified 
magnetic field close to the subshock is SB ~ 50/xG). This interpretation of the modest increase in the FWHM for 
low values of p max is further confirmed by the fact that the width is very weakly dependent on both £j n j and ?7th for 
Pmax S 10 TeV/c. 

For completeness, Figure [TU] shows the width of the intermediate component of the Balmer line as a function of the 
CR acceleration efficiency ecR, for V s h = 4, 000 km s -1 , ?7th = 0.5 and neutral fraction hiy = 0.5. The two curves refer 
again to /?down = 0.01 (upper curve) and /3d OW n = 1 (lower curve). In the upstream, energy is assumed to be deposited 
only in the protons, which therefore become warmer than the electrons. We assume that they do not equilibrate very 
effectively since the Coulomb collision timescale is typically longer than the advection one. 

The reason why the width of the intermediate Balmer line depends on the CR acceleration efficiency is that when 
€cr increases, the amount of TH close to the subshock increases as well, leading to a broader distribution of neutrals 
in the region affected by the NRF. The dependence of the width of the intermediate line on /?down lies in the fact that 
the NRF carries upstream the information about the ion temperature behind the shock. We further notice that the 
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Fig. 10. — FWHM of the intermediate Balmer line as a function of the CR acceleration efficiency, ecRi f° r a shock speed V s h = 4, 000 kms 
?7th = 0.5 and p max = 50 TeV. The solid and dashed lines refer to the cases ftjown = 1 an d /^down = 0.01 respectively. 



width of the intermediate component may also depend on the electron-ion equilibration upstream, (3 up . In Figure [TU] 
we assume no equilibration upstream, namely while ions are heated by NRF and TH, electrons always have T = 10 4 
K. Hence one has to be careful in using the intermediate line to infer shock properties, because this line is a function 
of several parameters. 

4. SUMMARY AND CONCLUSIONS 

In this paper we have presented the first formulation of the non-linear theory of diffusive shock acceleration at SNR 
shocks propagating in a partially ionized medium. At the same time, we have put forward the first fully self-consistent 
calculation of how the shape of Balmer H Q line is affected by the presence of accelerated particles. This development is 
crucial if one wants to use the width of the various components of the Balmer line to infer information on the efficiency 
of CR acceleration in SNRs. 

The three most important physical points discussed here in a quantitative way are the following: 1) neutrals pro- 
duced in CE events with hot ions downstream can return upstream and deposit energy and momentum there. This 
phenomenon, that we refer to as the neutral return flux (NRF), has an important dynamical effect on the structure 
of collisionless shocks as well as on the spectrum of accelerated particles (see also Paper I). 2) The structure of a 
collisionless shock is affected in a non-trivial way by the combined action of neutrals and CRs. 3) Also the shape of 
the Balmer H a line reflects the effects of the two phenomena mentioned above. 

The NRF was not predicted in previous calculations of the structure of SNR shocks, mainly because of their intrinsic 
limitation to fluid approaches or restriction to the upstream plasma alone. Both CRs and the NRF produce a precursor 
upstream of the shock: the CR- induced precursor develops on a scale of the diffusion length of particles at the maximum 
energy ~ D(p max )/V s h, while the neutral-induced precursor develops on a scale of a few CE interaction lengths. Both 
precursors cause heating of the ion plasma upstream. The heating associated to the NRF is mediated by CE and 
ionization upstream. The CR-induced heating is most effective only in the presence of turbulent (non adiabatic) 
heating. 

For typical values of the parameters relevant for SNRs, the spatial extent of the neutral-induced precursor is such as 
to overlap with the diffusion range of particles with p < 5—10 TeV/c. In this energy region the spectrum of accelerated 
particles is steepened by the dynamical effect of the NRF, as already found in Paper I, in the c ontext of test-particle 
calculations. The calcula tions presented here are a generalization of previous work on NLDSA (lAmato fc Blasil I2005L 



2006 ; iCaprioli et al.|[2009l ) and on the dynamical effects of neutrals on the structure of a collisionless shocks ( Blasi et al] 



2012t) . The technique adopted for the solution of the non-linear equations that describe the whole problem is basically 
an iterative method in two stages: in a first stage, we calculated iteratively the distribution of neutrals and ions 
(see Sections 12. 1[ 12.41 and Paper I). The distribution function at a generic location is not Maxwellian therefore the 
calculation of the neutral distribution is kinetic, while ions (protons and electrons) are treated as fluids. In a second 
stage, the updated profile of neutrals and ions is used to calculate the spatial distribution and spectrum of accelerated 
particles and waves (Sections 12.21 and |2.3|) . The procedure is repeated until convergence is reached. 

On average, as expected, the presence of neutrals leads to less concave and steeper spectra of accelerated particles 
than in a completely ionized medium. These effects are quite substantial for shocks with V s h < 3, 000 km s _1 . For faster 
shocks, the NRF becomes dynamically less important because neutrals are ionized downstream before they suffer a 
CE reaction. In these cases the shock modification is mainly induced by CRs. 

Besides the general interest in a theory of NLDSA in the presence of neutrals, we think that the present work is 
especially important in terms of providing an essential tool to interpret the shape of the H Q line and translate it into 
quantitative information on the presence of CRs in SNR shocks. The widths of the narrow and broad component of 
the Balmer line in SNR shocks have lon g been known to bear information on the temperature of the ions upstream and 
downstream of the shock respectively (|Chevalier fc Raymond! fl978t Ivan Adelsberg et all 11978ft . Both quantities are 
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sensitive, in turn, to the CR acceleration efficiency: if CRs arc efficiently accelerated, the plasma upstream is heated 
by the CR-induced precursor and the downstream plasma is heated less than in the absence of CRs because part of 
the ram pressure is channelled into CR acceleration. In terms of Balmcr line, one expects the narrow component to 
be somewhat broader and the broad component to be somewhat narrow er whenever CR accele ration is efficient. A 
phenomenological approach to the problem was recently put forward by IRavmond et al.l (|2011D . where however the 
CR spectrum and spatial distribution are both assumed rather than calculated. As described in the p resent paper, 
the CR properties are profoundly affected by the presence of neutrals. Moreover IRavmond et al.1 (|201lD do not treat 
the downstream plasma, therefore there is no NRF in their calculation. 

In Fig.[6]we can see the temperature of the plasma calculated self-consistently. Both the heating in the precursor (for 
different values of the efficiency of TH) and the reduced temperature in the downstream (as a result of CR acceleration) 
are clearly visible. 

These effects reflect in the shape of the Balmer line, which we calculated by using the formalism introduced in Paper 
II. The width of the narrow Balmer line is found to depend rather sensibly upon the level of TH and the value of 
maximum momentum p max . The dependence on p max can be understood, because increasing p max leads to larger CR 
precursors and therefore neutrals incoming from upstream infinity have a longer time to experience CE and adapt to 
the ion temperature. Larger TH reflects in an increased width of the narrow Balmer line. 

The level of TH does not affect appreciably the plasma temperature downstream, therefore the width of the broad 
Balmer line mainly reflects the CR acceleration efficiency. However, the broad line is sensitive to the level of equi- 
libration between protons and electrons downstream: if electrons and protons reach thermal equilibrium in a time 
short compared with CE (clearly this cannot happen through collisions) then the ion temperature is lower, thereby 
mimicking the presence of CRs. The estimate of the CR acceleration efficiency ecR might be less dependent upon this 
ambiguity if the width of the broad Balmer line is quite small. These cases however, suggesting a v ery high acceleration 
efficie ncy, are very difficult to obtain in the presence of TH and magnetic reaction on the shock (jCaprioli et al.l [20081 
120091 ). Moreover these cases would correspond to very concave spectra of accelerated particles and of the associated 
gamma rays, that so far have never been observed. The combined observation of the narrow and broad Balmer line 
might also help breaking the degeneracy between CR acceleration efficiency and clcctron-proton equilibration. By the 
same token, observation of an X-ray spectrum of thermal origin may independently provide a measurement of the 
electron temperature. 

Finally, we stress that the existence of the NRF leads to the formation of an intermediate component of the Balmcr 
line, with a width that depends on the temperature of protons downstream, on the electron-proton equilibration and 
on the level of TH. The practical possibility to infer ccr from data on individual SNR will need to be assessed on the 
case- by-case basis. 
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